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ABSTRACT 

A comparison of axisymmetric Direct Simulation Monte Carlo (DSMC) Analysis Code (DAC) 
results to analytic and Computational Fluid Dynamics (CFD) solutions in the near continuum 
regime and to 3D DAC solutions in the rarefied regime for expansion plumes into a vacuum is 
performed to investigate the validity of the newest DAC axisymmetric implementation. This new 
implementation, based on the standard DSMC axisymmetric approach where the representative 
molecules are allowed to move in all three dimensions but are rotated back to the plane of 
symmetry by the end of the move step, has been fully integrated into the 3D-based DAC code 
and therefore retains all of DAC’s features, such as being able to compute flow over complex 
geometries and to model chemistry. Axisymmetric DAC results for a spherically symmetric 
isentropic expansion are in very good agreement with a source flow analytic solution in the 
continuum regime and show departure from equilibrium downstream of the estimated breakdown 
location. Axisymmetric density contours also compare favorably against CFD results for the Ftl E 
thruster while temperature contours depart from equilibrium very rapidly away from the estimated 
breakdown surface. Finally, axisymmetric and 3D DAC results are in very good agreement over 
the entire plume region and, as expected, this new axisymmetric implementation shows a 
significant reduction in computer resources required to achieve accurate simulations for this 
problem over the 3D simulations. 


INTRODUCTION 

At high altitudes, plume impingement problems cannot be modeled using Computational 
Fluid Dynamics (CFD) solvers because the continuum assumption becomes invalid as the flow 
expands into a near vacuum. The Direct Simulation Monte Carlo (DSMC) method 1 , which models 
gas at the molecular level, is the method of choice to simulate such flows. Even at the low 
densities encountered in high altitude plume impingement problems, the actual number of 
molecules in the gas remains very large and cannot be simulated directly. Instead of simulating 
the actual gas molecules, the DSMC method uses representative molecules, usually between 
several millions to hundreds of millions, to model the gas flow. Representative molecules are 
tracked in space and time inside a meshed domain and the DSMC method assumes that the 
transport and collision phases can be modeled in two distinct computational steps. A network of 
computational cells is used to group neighboring molecules with possible collision partners and 
individual or groups of computational cells are also used to statistically sample the molecular data 
to produce macroscopic properties such a density and temperature. 

The main requirements to obtain an accurate solution using the DSMC method are 
threefold. First, the timestep size must be smaller than the local mean collision time to enable the 
decoupling of the transport and collision phases. Next the collision cell size must be smaller than 
the local mean free path, or mean distance travelled by a molecule between two successive 
collisions. This requirement ensures that artificial viscosity isn’t introduced into the solution by 
colliding molecules separated by more than one mean free path. Finally, each computational cell 
should contain a large enough number of representative molecules per cell to ensure an accurate 
statistical representation of the gas during the collision phase. 
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The DSMC method is well suited to solve plume expansions into a vacuum as it can 
theoretically model flows that span continuum, transitional and free molecular regimes. 
Unfortunately, based on the fundamental constraints of the method, accurate computations of the 
continuum regions of the flow with the DSMC method are prohibitively expensive. Additionally, 
the far field regions of the plume that are in the rarefied regime cannot be modeled using a 
continuum approach. Therefore, plume expansion problems are usually solved using a sequential 
approach. The near-field region of the expansion plume is usually computed using either a source 
flow model or using a CFD solver. Then, at a prescribed interface, the necessary data is 
extracted from the continuum solution and is provided as input to the DSMC code. Even when 
only modeling the far field region of the plumes the computational cost of the DSMC simulation 
may still be prohibitive, therefore, cost savings approaches, such as taking advantage of the 
symmetry of the problem, whenever possible, must be considered. Examples of plume 
impingement problems solved using the DAC code include Space Shuttle’s airlock venting gas on 
the Hubble Space Telescope (HST) solar array 2 , high fidelity plume impingement on the 
International Space Station (ISS) 3 , and Orion self-impinging plumes during re-entry 4 . 

The purpose of the present analysis is to validate the axisymmetric approach that has 
recently been implemented in NASA’s DSMC Analysis Code (DAC) for the simulation of high 
altitude plumes. Axisymmetric results are compared to an analytic solution for a spherically 
symmetric gas expansion into a vacuum and to CFD and to 3D DAC results for the R1 E thruster 
plume. 


SIMULATION METHOD 

One of the primary objectives in developing NASA’s DAC code 5 was to provide a high 
fidelity modeling tool for vacuum plume impingement on spacecrafts and satellites. 

PLUME SIMULATIONS 


Simulating plume expansions into a vacuum can be a challenging problem as no single 
computational approach can accurately and/or efficiently model the gas from the combustion 
chamber to the far field region where the gas impinges on a given surface. Instead, a multi step 
approach is generally chosen where different methods appropriate for different regimes of the gas 
flow are used sequentially. For a given thruster, a CFD approach is used to model the flow inside 
the nozzle and into the near field region of the thruster. At some distance from the nozzle exit 
plane, the gas is sufficiently expanded that the continuum assumption breaks down and the CFD 
solution becomes invalid. The limit between continuum and non-continuum regions of the flow is 
usually described using a Bird breakdown parameter, B, equal to 0.05, where B is given by 6 : 


where A is the mean free path and c is the average molecular speed of the gas. 

The computed interface is generally smoothed out and is connected to the nozzle exit lip 
to form the inflow surface that will be used in the DSMC computation. The last step of this 
process is to interpolate the CFD results to the inflow surface such that required input, such as 
density, temperature and gas velocity, can be provided to the DSMC code. 


DAC AXISYMMETRIC IMPLEMENTATION 


The DAC code is a fully three dimensional parallel DSMC code that has been used over 
the years to solve a wide range of rarefied problems from high altitude re-entry flows to plume 



impingement studies. In order to solve flows around complex geometries, the surfaces are 
represented in DAC by triangulated, topologically water tight meshes that are then embedded 
inside a two-level Cartesian grid system. The grid system can be refined following the DSMC 
requirements on cell size and timestep size using an iterative adaptation process. The DAC code 
models the relevant energy exchange between the translation, vibrational and chemical modes of 
the molecules and has been designed to support fully three dimensional, two dimensional, and 
axisymmetric simulations. 

While inherently three-dimensional, due to the random movement of the molecules, the 
standard DSMC method can be modified to solve axisymmetric problems. The most common 
axisymmetric implementation allows molecules to move in all three dimensions but forces them 
back to the symmetry plane by the end of the transport phase. One drawback of this approach is 
that, depending on the actual implementation, this approach may produce depleted density and 
pressure contours near the axis of symmetry due to the fact that the cell volumes are very small 
in that region of the domain 7 . The original DAC axisymmetric implementation used a different 
approach that employs a physical wedge to constrain molecules to a small region of space 
around the symmetry plane. Unfortunately, this approach could not guarantee that the cell size 
was smaller than the mean free path everywhere in the domain, more specifically in the regions 
far away from the axis of symmetry. The new DAC axisymmetric implementation which is a hybrid 
between the original method used in DAC and the most common method is presented in more 
detail in the present section. 

Because DAC was initially implemented as a 3D code, the present approach requires the 
input of a 3D surface geometry (Figure 1). For a plume simulation, the outline in the symmetry 
plane of the breakdown surface for the CFD solution is extruded in the third direction to create a 
triangulated inflow surface where molecules will be input to the DSMC simulation. This surface is 
integrated with other walls, modeled as outflow surfaces for plume simulations, that will help 
define the computational domain as shown in Figure 1 . As long as the extracted surfaces extend 
outside of the domain limits for the Cartesian mesh, the surface geometry is no longer required to 
be water tight. This approach provides a way to save on memory usage as fewer triangles are 
required to define the surface geometry. 


Outflow walls 



Figure 1 : Schematic of the surface geometry for a DAC axisymmetric simulation of a 
spherically symmetric expansion plume into a vacuum. 


Because the surface geometry and the Cartesian grid no longer model the actual 
symmetry of the problem, a pre-processing step has been added to the DAC code to accurately 
model the axisymmetric problems. First, the current approach computes the cell volumes in three 
dimensions and prost-processes them to obtain the actual axisymmetric cell volumes. For both 
computational cells with and without intersecting surfaces, the present approach computes the 
axisymmetric volumes to within a few percent of the actual values. Next, the surface triangles are 
intersected with the symmetry plane in order to produce inflow and solid line segments. Because 
of the symmetry of the problem, each line segment may represent a different rotated shape and 
this has to be taken into account during both the creation and the move phases. The vertical 
segments represent flat rings, the horizontal segments represent cylinders, while all other 
segments represent conical sections. Molecules are created from the inflow line segments and 
enter the domain inside the symmetry plane. Next, the molecules are moved in all three 
dimensions. Because DAC allows variations in timestep and ratio of real to simulated molecules 
from one Level-1 cell to another, the code first computes the time it takes a molecule to reach the 
edge of the cell. The main difference from the 3D solver is that in the direction perpendicular to 
the axis of symmetry, the code must compute the time it takes the molecule to reach the cell edge 
after it has been rotated back to the symmetry plane. If the molecule crosses to another cell, the 
timestep size required to reach the edge is computed and the molecule is moved by that 
timestep. Next, the code checks if the 3D molecule path intersects any of the surfaces 
represented by the solid line segments and, if it does, the code computes the time it takes the 
molecule to reach the surface and models the surface interaction. At the end of this sub-step, the 
molecule is rotated back to the symmetry plane. This process is repeated until the molecule has 
been moved by the full timestep. The other steps in the DSMC method, such as the molecular 
collision and sampling phases, remain similar to the 3D implementation. Once the simulation is 
over, surface properties can be computed and are output in line plot format after being possibly 
smoothed out or interpolated to uniform segments. 


RESULTS 

SPHERICALLY SYMMETRIC ISENTROPIC EXPANSION INTO A VACUUM 


A spherically symmetric flow that isentropically expands into a vacuum is an example of a 
self-similar flow because the gas properties only depend on the radial distance from the point of 
origin. Using the isentropic relationships and number density, n*, and temperature, T*, at the sonic 
radius, r*, the gas properties at radius r are given by: 
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where y is the ratio of specific heat, R is the ratio of the ideal gas constant by the gas molecular 
weight, and M, T, V, and n are the Mach number, temperature, velocity and number density at 
radius r, respectively. 




Unfortunately, this analytic solution is only valid in the near field region of the expansion 
because, as the gas expands, the flow becomes rarefied and the continuum assumption used to 
produce Equation (2) breaks down. The objective of the present simulations is to show that the 
DAC axisymmetric implementation correctly models the physics of the expansion in the 
continuum regime while also accurately modeling the departure from equilibrium as the flow 
becomes rarefied. Therefore, the interface at which the analytic solution is passed as input to the 
DAC simulation is chosen to be in the continuum region of the flow and instead of specifying the 
number density at the sonic radius, the Bird breakdown parameter is used. 

For the present self-similar flow, the Bird breakdown parameter at radius r can be 
rewritten as: 


B = 2 
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Replacing the mean molecular speed and using the hard sphere approximation for the mean free 
path, the number density can be expressed as function of the Bird breakdown parameter as: 
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where d is the molecular diameter. 


The present simulation models the expansion of an argon perfect gas with a ratio of 
specific heats of 5/3, a molecular weight of 39.948 g/mol and a collision diameter, d, of 4.17e 10 
m. The sonic conditions occur at a radius of 1 m from the point of origin of the expansion with a 
gas temperature of 1000 K and a Bird breakdown parameter of 0.0025. The inflow surface where 
molecules are created in the DAC domain is located at a radius of 1 .05 m. As shown in Figure 2, 
the flow conditions at the inflow surface are still in the continuum regime and the continuum 
assumptions are expected to break down at a radius of about 2.5 m. The DAC simulation is 
started on a uniform grid with a cell size of the order of the mean free path for a gas at a number 
density one order of magnitude smaller than the gas at the inflow surface. Because this grid 
under-resolves the flow near the inflow surface, the simulation uses virtual sub-cells in order to 
reduce the effective collision cell size 8 . The total number of molecules in the domain is 
approximately 1 75 million and cells near the inflow surface have as many as 1 0,000 molecules. 
This simulation took less than one hour to run using 120 processors. Based on this initial 
simulation, an adapted grid is created with a grid size everywhere in the domain of the order of 
the local mean free path. Because the number density in the entire domain is actually much lower 
than the value used to create the uniform grid the adapted simulation uses a lot fewer cells and 
therefore molecules, only approximately 4 million, while still following the basic DSMC constraints 
on cell size, timestep size and number of molecules per cell. As a consequence, the adapted 
case is much faster than the uniform case and takes only five minutes to complete when using 
120 processors. 

Axisymmetric implementations in DSMC codes have been shown in the past to have 
problems simulating the gas flow near the axis of symmetry because the cell sizes in that region 
of the domain become very small 7 . The present axisymmetric number density contours (Figure 2), 
however, show no depletion near the axis of symmetry and the contours are concentric over the 
entire domain. In order to directly compare the DAC axisymmetric and the analytic solutions, data 
have been extracted along the white lines shown in Figure 2. The three different lines, in the +X 
direction, in the +Y direction and along the Y = -X direction, were chosen in order to validate the 
solution over the entire domain. The non-dimensional number density DAC line plots as function 
of radial distance from the point of origin (Figure 3a) are in very good agreement over the entire 
domain with the analytic solution. The non-dimensional temperature line plots, however, show 
that for radial distances greater than around 5 m, the DAC temperatures depart from the 
equilibrium temperature obtained with the analytic solution. As the gas becomes more and more 
rarefied, the local collision rate decreases to the point where there is no longer enough collisions 
for the gas to remain in local thermal equilibrium. The temperature freezes and therefore departs 
further and further from the analytic temperature which continues to decrease with radial 
distance. This departure from equilibrium happens downstream of the location at which the Bird 



breakdown parameter is equal to 0.05 which means that for this case the initial assumption that 
the continuum assumption breaks down at that location is conservative. 
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Figure 2: DAC axisymmetric number density contours for a spherically symmetric 

expansion into a vacuum. 
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Figure 3: DAC and analytic non-dimensional number density and temperature line plots as 

function of radial distance. 



R1E THRUSTER PLUME 


A more realistic test of the new axisymmetric implementation is to simulate an actual 
thruster plume expansion into a vacuum. The R1 E thruster is used in the present simulations to 
compare the axisymmetric DAC results to both 3D DAC results and a CFD solution. The nozzle 
section of the thruster and the near field region of the plume have been computed using a CFD 
solver. Figure 4 shows the Bird breakdown parameter contours for the CFD solution. Please note 
that the contours are not matching at the boundaries between the CFD sub-domains because 
gradients aren’t computed based on data located on both sides of the boundaries. The inflow 
surfaces for the axisymmetric and 3D solutions are created based on the contour line where the 
Bird breakdown parameter is equal to 0.05 and by respectively extruding the line in the Z 
direction or rotating the line 360 degrees. For reference, the outline of the inflow surface in the 
symmetry plane used in the DAC simulations is also shown in Figure 4. Once the inflow surface 
has been created, the CFD solution is interpolated to each surface. In the present DAC 
simulations, the thruster gas is assumed to be a single species with a molecular weight equal to 
the average molecular weight in the CFD solution at the interface near the centerline of the 
plume. The axisymmetric simulation ran under thirty minutes on 48 processors and an iterative 
adaptation process was again used in order to properly refine the mesh and timestep size 
everywhere inside the computational domain. For the adapted grid simulation, the grid resolved 
the local mean free path everywhere in the domain and the number of molecules decreased to 7 
million from the 8 million used in the uniform grid simulation. The 3D run took about nine hours to 
run on 1 80 processors and required 90 million molecules for the simulation on the uniform grid 
and 630 million molecules for the simulation on the adapted grid. The simulation on the adapted 
grid was also resolved to mean free path resolution by using virtual sub-cells for selecting the 
collision pairs 8 . 
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Figure 4: CFD Bird parameter contours in the symmetry plane for the R1E thruster. 


Non-dimensional density and temperature contours are shown in Figure 5. Overall, the 
3D and axisymmetric DAC density and temperature contours are in very good agreement over 
the entire domain. The DAC density contours also compare favorably to the CFD solution except 
very near the left boundary of the domain. In the CFD solver, the region near the nozzle lip is 
simulated as a viscous wall. In the DAC simulation, however, there is no physical wall there and 
the gas is allowed to freely expand. This difference in the modeling of the boundary can be seen 



in the density increase near the wall for the CFD solution. As with the expansion, when 
comparing the DAC solutions to the CFD results, non-equilibrium effects are mostly apparent in 
the temperature contours. The non-dimensional DAC temperature contours are very different 
from the CFD contours except very near the axis of symmetry, near the inflow surface. Therefore, 
for this case, the departure from equilibrium occurs near the location at which the Bird breakdown 
parameter is equal to 0.05. Similar to the expansion flow, the temperature computed in the DSMC 
simulations freezes as the gas expands and is therefore higher in most of the domain than the 
temperature computed in the CFD simulation. 
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Figure 5: Non-dimensional density (left) and temperature (right) contours in the symmetry 

plane for the R1 E thruster. 


SUMMARY AND CONCLUSIONS 

Axisymmetric DAC results have been compared to an analytic solution for a simple 
spherically symmetric expansion. The results are in very good agreement with the analytic 
solution in the continuum region of the domain and show the expected departure from equilibrium 
in the temperature line plots downstream of the estimated continuum breakdown location. 
Axisymmetric results have also been compared to CFD and 3D DAC solutions for the R1 E 
thruster. For this case, the density contours are in very good agreement with the CFD contours 
except very near the left boundary of the domain and freezing of the temperature is observed 
everywhere in the DSMC domain. Finally, the axisymmetric and 3D DAC solutions are in very 
good agreement over the entire domain and as expected the axisymmetric results were 
generated faster than the 3D results by more than one order of magnitude. 
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